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Abstract 

The solvation force of a simple fluid confined between identical planar walls is studied in two 
model systems with short ranged fluid-fluid interactions and long ranged wall-fluid potentials de- 
caying as — Az~ p , z — ► oo, for various values of p. Results for the Ising spins system are ob- 
tained in two dimensions at vanishing bulk magnetic field h = by means of the density-matrix 
renormalization-group method; results for the truncated Lennard-Jones (LJ) fluid are obtained 
within the nonlocal density functional theory. At low temperatures the solvation force f so i v for the 
Ising film is repulsive and decays for large wall separations L in the same fashion as the boundary 
field f so i v ~ L~' p , whereas for temperatures larger than the bulk critical temperature f so iv is attrac- 
tive and the asymptotic decay is f so i v ~ L~( p+1 \ For the LJ fluid system f so i v is always repulsive 
away from the critical region and decays for large L with the the same power law as the wall-fluid 
potential. We discuss the influence of the critical Casimir effect and of capillary condensation on 
the behaviour of the solvation force. 

PACS numbers: 
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I. INTRODUCTION 



Fluid mediated interactions between two surfaces or 

n 



arge particles, usually referred to 
2 1 that can be relevant in colloidal 



nonolavers 



as solvation forces f so i v may lead to subtle effects 
systems and for many applications and modern technologies such as lubrication, adhesion 
or friction. There is a rapidly growing literature on experimental investigations and other 
phenomena associated with solvation forces. Direct measurements with the Surface Force 
Apparatus (SFA) 3J have revealed the richness of behaviour of these forces, for example, 
their sensitivity to specific properties of the intervening fluid. Dependence of the solvation 
force on the chemical and physical properties of confining surfaces is also of much interest. 
Different surfaces have been developed and used in SFA by adsorbing or depositing a thin 
film of some other material on a mica surface, for example, lipid monolayers or bilayers, 
metal films, polymer films or other macromolecules such as proteins 

On the other hand, the theory of the solvation forces which would help to interpret mea- 
surements performed for different liquids and different surfaces is far from being complete 
even for simple model systems. One of the issues which has not been systematically in- 
vestigated is how the properties of the solvation force f so iv(L), such as its sign, magnitude 
and asymptotic dependence on the distance between surfaces, vary with the range and the 
strength of the substrate-fluid potential and with the thermodynamic state. In the present 
paper we address these issues for a simple fluid confined between two identical parallel solid 
substrates separated by a distance L. 

The majority of results available for the solvation force in this simple system comes from 
a model fluid in which both the fluid-fluid interatomic potential and the wall-fluid potential 
V s (r) are of short range. Here, by f so iv{L) we are referring only to the fluid contribution to 
the force, the direct wall-wall contribution is not included. Let us summarize briefly results 
that are the most relevant for the present paper. The precise form of the solvation force 
depends on the details of V s (r), the bulk state point, as well as on L. Theory and computer 



simulations show that at small wall separations the solvation force exhibits oscillatory be- 
haviour with L. These oscillations, observed in direct measurements of the solvation force in 
real fluids, reflect packing effects that produce oscillations in the density profile for liquids 
near walls 7|. The behaviour of f so iv{L) reflects the behaviour of the density profile also 
in the limit of large separation, L —>■ oo. General statistical mechanical arguments predict 
that far from the bulk critical temperature T c and from any phase transitions (condensation 
phenomena) f so i v {L) for L — > oo must decay to zero in the same fashion as the profile at a 
single wall decays to its bulk value, i.e., as the radial distribution function g{r) of the bulk 



fluid |S, 



9 



1(1 lllj. Thus, for V s (r) of finite range and L — > oo, f so i v {L) ~ e~ a ° L or, for suf- 
ficiently high densities and/or temperatures, (states on the oscillatory side of the so-called 
Fisher- Widom line) f so iv(L) ~ exp(—a L) cos(aiL) . The quantities a and a.\ depend on 
only the bulk pair correlation function. The sign of the solvation force and its temperature 
dependence at fixed L are available from mean-field (lattice, Landau, density functional) 



analyses 
strips 
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14, 
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and from exact transfer matrix methods for two-dimensional 
19j . For simple fluids confined between identical walls f so iv(L) < for large L, i.e., 
the net force between the substrates is attractive for large separations. Away from the bulk 
critical temperature T c , f so i v (T) is small in magnitude (weakly attractive). On approach- 
ing T c at vanishing ordering field, h — 0, where h ~ [A — fi sa t, where fi sat is the chemical 
potential at coexistence, the force increases in magnitude and exhibits a minimum. For 
strongly adsorbing walls this minimum is located above T c and occurs for L of the order 



of the bulk correlation length Critical scaling arguments predict 
criticality f so iv(L) decays algebraically for large separations, i.e., 
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2l| that at bulk 



fsoiv{L) ~k B T c A 12 (d-l)L- 



oo. 



(1) 



d is the spatial dimension and a universal number (Casimir amplitude) for 2 < d < d>, 
where the upper critical dimension d > = 4 for the Ising universality class. A12 is negative 
for identical walls, i.e. the Casimir force is attractive. 



In the present paper we study two model systems in which fluid-fluid interactions are 
short ranged but the substrate-fluid potential V s (r) is long ranged, and enquire how the 
features of the solvation force described above change with the range and the strength of 
V s (r). As a first system we consider a (Ising) lattice gas in a slit geometry subject to identical 
boundary fields which depend on the distance I of a lattice site from the boundary 



VJr) = Hf 



hi 
IP 



(2) 



with h\,p > 0. The Ising model is particularly useful in fundamental studies of finite-size 
and surface effects in confined fluids since exact calculations are available, at least in two 
dimensions and for short ranged boundary fields 
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24, 
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. Another advantage 



is that in d = 2 this model is amenable to the systematic investigation for arbitrary boundary 



and bulk fields by means of density- matrix renormalization-group (DMRG) method [2jj . The 
DMRG method, based on the transfer matrix approach, provides a very efficient algorithm 
for constructing the effective transfer matrices for large systems. Comparisons with exact 
results for the case of vanishing bulk magnetic field and boundary fields acting only on 
spins in the surface layers show that this technique gives very accurate results in a wide 
range of temperatures, including near the bulk critical temperature, and for large widths 



of the strip 
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29j. Here, using the DMRG method we calculate f so i v as a function 
of temperature along the bulk two- phase coexistence line, i.e., for vanishing bulk magnetic 
field, for different values of p, hi and L. We also study the asymptotic decay, with L, of the 
solvation force for strongly adsorbing walls (large hi) and temperatures away from T c . 

Lattice models of a confined fluid cannot describe packing effects that are reflected in 
the oscillations of the solvation force at small wall separations. The other failing of the 
(Ising) lattice gas model is that it has an exact particle-hole symmetry. For real fluids such 
symmetry is only approximate which can be of relevance for the properties of f so i v away 
from criticality. Therefore, to make our analysis more complete we also study a Lennard- 



Jones (LJ) fluid in a slit geometry. Specifically, we consider a truncated LJ 12-6 fluid-fluid 



intermolecular pair potential and a wall-fluid potential of the form: 



V s (r) = V s (z)=Ae fw 



2 /ct\ 9 /aV 
15 \z) ~ \z 



(3) 



where p = 2,3; a is the diameter of the fluid species, while Sf w describes the strength of 
the wall-fluid interactions. We note that for p = 3, V s (z) models a wall which is assumed 
to comprise a half space of LJ particles 4. We calculate f so i v using a nonlocal density 
functional theory (DFT) along a similar thermodynamic path as in the Ising system, i.e., as 
a function of the temperature along the bulk two phase-coexistence line and at the critical 
density for T > T c . Our results refer to the densities on the liquid side of this line. We also 
investigate the asymptotic dependence of the solvation force on the distance between the 
walls at fixed temperature away from T c . 

The asymptotic behaviour of the solvation force for the Lennard- Jones fluid follows from 



the analysis by Attard et al. 3J| based on the wall-particle Ornstein-Zernike (OZ) equations. 
Using the hypernetted-chain closure these authors derived the interaction free energy per 
unit area between the planar walls F 0Q (L) as a convolution of wall-solvent pair-correlation 
functions. For a power-law fluid-fluid interaction potential — Ar~ n ,n > 3,r — > oo and a 
wall-fluid potential decaying as —Bz~ p for z —>■ oo a formula for the asymptotic behaviour 
ofF 00 (L) (Eq.(6.11) ofRef.Q) is 

Foo(L) ~ (3u 00 (L) + ^-L^ - —L A ~ n , L - oo (4) 

p — 1 [n — 2) [n — 3) [n — 4) 

where Uqo(L) is the direct wall-wall interaction potential per unit area and p is the density of 
the fluid. From Eq. (0J) it follows that in the case of truncated LJ fluid, when the power-law 
tails are omitted, the behaviour of f so iv(L) = — (dF 00 (L)/dL) T for large L is 

fsoiv(L) ~ ^j^-, L -» oo, (5) 

where B — Aef w and we have neglected the contribution due to the direct wall- wall interac- 
tion potential. Thus the solvation force is repulsive and decays with the same power law as 



the wall-fluid potential. The above prediction is treated as a reference point for the analysis 
of the asymptotic behaviour of our results. We find that, as one expects, Eq. (JSJ) is indeed 
valid provided one is away from the critical temperature and from any phase transitions. 
The influence of the critical Casimir effect and of capillary condensation on the behaviour 
of the solvation force is also discussed in the present paper. 

It is reasonable to expect that away from T c and from any phase transitions the 'magnetic' 
analog of the solvation force in system of Ising spins subject to identical boundary fields 
should have the same asymptotic decay power law as the boundary field Hf. Somewhat 
surprisingly, this is not the case for temperatures far above T c where we find that f ao iv(L) 
decays as for L — > oo. For low temperatures the decay agrees with the prediction ©. 

In order to understand the nature of this particular behaviour we analyse the appropriate 
Landau theory. 

The paper is organized as follows. Section HH is devoted to the Ising model studies. 
We define the model and briefly review the known results pertinent to the present studies. 
We then proceed to describe the results of the DMRG studies. In Sec. IIHI a continuum 
Landau theory is investigated for long ranged surface fields. In Sec. |TV] DFT calculations 
are presented and discussed. Section IVl summarizes our results and makes some conclusions. 

II. ISING MODEL RESULTS 
A. The model 

We consider an Ising ferromagnet in a slit geometry subject to identical boundary fields. 
Our DMRG results refer to the d = 2 strip defined on the square lattice of the size L x 
M, M — > oo. The lattice consists of L parallel rows at spacing a = 1, so that the width of 
the strip is La = L. We label successive rows by the index I. At each site, labeled (l,k), 
there is an Ising spin variable taking the value o~ik = ±1. We assume nearest-neighbour 
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interactions of strength J and Hamiltonian of the form 

L 

H = -J{ vikVi'k' + ^HiJ^vik}- (6) 

<lk,l'k'> 1=1 k 

The first term in © is a sum over all nearest-neighbor pairs of sites while in the second 
term Hi = Hf + if£ +1 _j is the total boundary magnetic field experienced by a spin in row 
I. The single-boundary field iff is assumed to have a form: 

with p > 0, and the reduced amplitude of the boundary field hi > 0. 
B. Solvation force for short ranged boundary fields 

All previous work on the behaviour of the solvation force in Ising films used localized 
boundary fields acting only on spins in surface layers I = 1 and I = L: 

Hi = Mi,, + h 2 5 L>l . (8) 

In the limit hi,h,2 — > oo Eq. (jHJ corresponds to the (++) fixed spin boundary conditions, 
widely studied in the literature. 

The total excess free energy per unit area for the case of identical surface fields hi = h 2 
and vanishing bulk magnetic field h = can be written as 

f ex (L) = L(f(L, T, h) - f b (T)) + 2f w (T, hi) + f*(L, T, hi) (9) 

where / is the total free energy per site, /& is the bulk free energy, f w is the L-independent 
surface excess free energy contribution from each surface, and /* is the finite-size contribution 
to the free energy. All energies are measured in units of J and the temperature in units of 
J/ks- /*, which vanishes for L — > oo, gives rise to the generalized force, which is analogous 
to the solvation force between the walls in the case of confined fluids 

fsoiv = -(df ex (L)/dL) TM . (10) 



For identical surface fields the solvation force is attractive for all thermodynamic states, i.e., 
fsoiv < 0. The asymptotic decay of the solvation force depends on the temperature range. 
For temperatures sufficiently far away from the bulk critical temperature T c , f so i v decays as 



exp(— L/£b), where is the bulk correlation length 



as a resu. 
effect 



t of critical fluctuations 



1 



Near T c , f so i v becomes long ranged 



a phenomenon which is known as the critical Casimir 
At T = T c and h = the asymptotic decay is given by Eq. (JTJ). For the case 
of (++) boundary conditions the temperature dependence of the solvation force at fixed L 



and h = 0, or equivalently the scaling function W ++ (y) is defined by 



fsoiv/k B T c = (d-l)L- d W ++ (y), 



(11) 



where the scaling variable y = t{L/^q) 1 ^ . In d = 2 the scaling function was determined 
by exact transfer matrix methods [19]. Here t = (T — T c )/T c and £ = £o r_I/ (with v = 1 
in d = 2) is the bulk correlation length. It was found that f so i v plotted as a function of 
temperature attains a pronounced minimum above the bulk critical temperature T c when 
y = t(L/£q~) 1 / 1 ' = 2.23, or L ~ 2.23£. The amplitude of the correlation length above T c is 
£ + w 0.5673. 

It was also shown 



191 ] that the scaling function has the property, for d — 2, 



(12) 



where subscript 00 refers to hi = hi = 0. This implies that at h — and fixed L the 
function f so iv(T) evaluated for free boundaries has its minimum below T c . The behaviour of 
the function f so i v {T) in the crossover between hi = oo and h\ — has not been studied. In 
the next section we perform such an investigation for both short ranged and long ranged 
boundary fields. 

Finally, we note that for free boundaries h\ = hi = the location of the minimum of 
the solvation force, at h — 0, is associated with the critical temperature, T r r„ which for 



d > 3 and large but finite L denotes the end of the two-phase coexistence 



19 



32j. T c r 



lies on the h = axis and is shifted below the temperature of the bulk critical point T c 



by an amount given by the expression 33 1, r ~ —L~ x l v . For nonvanishing surface fields, 
hi = h>2 > 0, the situation is different. In this case the preferential adsorption of (+) spins 
at each wall leads to a shift of the bulk phase boundary in the (h,T) plane to h < 0. 
This phenomenon of capillary condensation strongly influences properties of the Ising films 



- also above the (capillary) critical point (h C: L(hi),T C;L (hi)) < T c 



32, 



341 ] . The minimum of 



the scaling function W ++ (y) lies above T c which can be accounted for by the fact that the 
most pronounced features in the solvation force occur along the continuation to higher T of 

n 

the capillary condensation line h co (T) [32]. The same is true for other thermodynamic and 
structural quantities that arise in the film geometry, such as the specific heat, adsorption, or 
longitudinal correlation length £n. Specifically, f so i v at fixed T > T c L has a deep minimum 
at some h < 0, which corresponds roughly to the continuation of the line h co (T). As the 
temperature increases the minimum approaches h = and decreases in depth. The stronger 
hi, the bigger is the shift of the capillary critical point from the bulk coexistence line and the 
further above T c reaches the line along which the most pronounced minima of the solvation 
force occur. The shape of the scaling function W ++ {y) reflects this behaviour. W ++ (y) is 
very weak for T < T c when the condensation line is far away from h = 0, and develops the 
minimum above T c when the minima of the solvation force that occurs along the continuation 
of h co (T) approach the line h = 0. 

C. DMRG results for long ranged boundary fields 

The density-matrix renormalization group (DMRG) method was introduced in 199 2 by 
White as a numerical algorithm to study ground-state properties of quantum-spin chains 35] . 
In spite of the name, the method has only some analogies with the traditional renormalization 
group being essentially the numerical, iterative basis, truncation method. Later, the DMRG 



was adapted by Nishino for two-dimensional classical systems at non-zero temperatures 21 

10 



It is particularly well suited to study systems in confined geometry since it deals naturally 
with lattices of a size Lxoo. Generally, the DMRG method works best with open boundary 
conditions, which makes the technique appropriate to take into consideration the effects of 
surfaces. We have implemented a finite-size version of the DMRG algorithm designed for 



accurate studies of finite-size systems 



For a comprehensive review of a background, 



achievements, and limitations of the method, see Ref. 



36] 



In a transfer matrix approach a leading eigenvalue Xl of a transfer matrix Ti 

T L \v L ) = \ L \v L ) } (13) 

gives a free energy per spin of an Ising strip as 

Pf{L) = ~\n\ L . (14) 

The components of the eigenvector \vl) related to the leading eigenvalue give probabilities 
of various configurations. For classical spin systems the DMRG method is based on the 
transfer-matrix approach, where the leading eigenvalue and its eigenvector of the effective 
transfer matrix are calculated numerically. This method can be employed for a number of 
problems for which no exact solutions are available (e.g. Ising systems in a presence of bulk 
magnetic field, or, as in the present case, subject to long ranged surface fields). 

The main idea of the DMRG technique is to avoid the proliferation of states when the 
size of the system grows. Generally, a number of configurations in a Hilbert space of an 
Ising strip grows very fast with its width L (as 2 L ). Therefore, it is practicaly impossible 
to solve exactly systems with L > 25. In the DMRG approach one eliminates the least 
probable (in the density matrix sense) states and keep only the most important ones. From 
this stage calculations are not exact anymore, but we can obtain a very efficient approach 
if the weights of the discarded states are very small. Starting with a small system (e.g. 
L = 4 in our case), for which Tl can be diagonalized exactly, one adds iteratively two spin 
rows in the middle of a strip until the allowed (in the computational sense) size of effective 

11 



matrices is reached. Then further addition of new spins forces one to discard simultaneously 
the least important states to keep the size of effective matrices fixed. This truncation is 
done through the construction of the reduced density matrix whose eigenstates provide an 
optimal basis set. The size of the effective matrix is then substantially smaller than the 



original dimensionality of the configurational space [37[ (2m) 2 <C 2 L . Generally, the larger 
is m, the better is the accuracy. In the present case we keep this parameter up to m = 50. It 
is worth mentioning that at low temperatures in order to renormalize a transfer matrix we 
have to use its two eigenvectors corresponding to phases with the opposite magnetization. 
To construct the reduced density matrix from the lowest eigenstates one has to diagonalize 
an effective transfer matrix Ti at each DMRG step. Therefore, we used the so-called Arnoldi 
method Q. 

To calculate the solvation force we proceed in the same way as in the case of the short 
ranged surface fields. We calculate the excess free energy per unit area L&^L) = (/ — /&) L 
at Lq + 2 and Lq. For vanishing bulk magnetic field /& is known exactly 39]. Having values 
f ex (L + 2) and f ex (L ) we approximate the derivative in Eq. (jlOj) by a finite difference 
f solv = -(l/2)(/ e *(L + 2) - f ex (L )). 

We have chosen four different power- laws describing the decay of the surface fields, i.e., 
we consider Eq. (J2J) with p = 50, 4, 3, 2. The behaviour of the system with p = 50 should be 
close to the one with the short ranged surface fields, p = 4 corresponds to dispersion forces 
in d = 2. For each type of boundary field we calculate the solvation force as a function of 
the temperature f so iv{T) along the bulk coexistence line h = for a range of amplitudes 
hi. Results for a fixed strip width L = 300, and fixed field strength parameter hi ~ 8.16 
are plotted as a function of the scaling variable y = rL x l v (see Eq. (111)1 ) in Fig. ^ For 
this value of hi the system with short ranged surface fields is almost identical to the infinite 
surface field scaling limit, i.e., to (++) boundary conditions. The inset shows a magnified 
plot of the region around the bulk critical temperature T c w 2.269185. We have found that, 



12 



as expected, for p = 50 the solvation force behaves as for the case of short ranged surface 
fields, i.e., it is always negative, exhibits a minimum located above T c for y m 2.23, and away 
from the critical region approaches zero from below. For all other values of p we observe 
that f so i v is repulsive below T c . As T approaches T c the critical fluctuations give rise to the 
Casimir effect. The fluctuation- induced Casimir force is always attractive and we can see 
that f so i v changes its sign to become negative as the temperature increases towards T c . This 
effect becomes stronger on increasing the range of the surface fields (see Fig. [Q. Strikingly, 
for high temperatures f so i v remains attractive and approaches zero from below independently 
of the value of p. Notice, that for p = 2 the solvation force is distinctly stronger than for 
the higher values of p. 

Turning now to the asymptotic behaviour of the solvation force we plot in Fig. El In f so i v 
versus ln(l/L) for various L between 20 and 320. The force has been calculated for fixed 
hi ~ 8.16 at three different temperatures, corresponding to T C T c , T = T c , and T ^> T c . 
For p =50 we find asymptotic decay typical of that for short ranged surface potentials, i.e., 
away from T c the solvation force decays exponentially with L. In Fig. El we display only 
results for the parameter p = 4 and 2 but for p =3 we observe the similar behaviour. Filled 
symbols in FigEt represent results calculated at T/T c = 0.79, circles for p = 2 and squares 
for p = 4. The straight lines in this figure have slopes equal to p (the solid line has a slope 2 
and the dashed line has a slope 4) and fit the data at this temperature, far below the critical 
temperature. Symbols in Fig|2b represent results calculated at T/T c = 1.23. The straight 
lines in this figure have slopes equal to p + 1 (the dotted line has a slope 3 and the dot-dashed 
line has a slope 5). Symbols not connected by lines in FigEH are results obtained at the bulk 
critical temperature. The critical scaling analysis predicts that if the long ranged boundary 
field decays sufficiently rapidly, i.e., when (1/2) {d + 2 — rj) — p is negative, the boundary 
field is irrelevant in the RG sense 40j. r] is a critical exponent of a spin-spin correlation 
function. The correction to the leading short-ranged behaviour given by (jll|) should be of 
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41 1 . In the present 



the order of £-(<*-2+»?)/2-p and may become dominant for L 3> ^log^ 
case of the d = 2 Ising model rj is equal to 0.25 so that for p > 2 we expect to see the power 
law L~ 2 for the asymptotic decay of the solvation (Casimir) force at T = T c (see Eq. (JT|)). 
We have found very good agreement with this prediction for all values of p. One can see 
in Fig. |2t that both circles and squares corresponding to p = 2 and 4, respectively, align 
almost parallel to a straight line of slope 2. We have checked that for p = 4 the data fit 
a line with a slope ~ 2.02, whereas for p = 2 the fitted line has a slope ~ 2.14. Notice 
that for p = 2 the correction to the leading decay of the solvation force is of the order of 
£-2.125 j n orc [ er to observe a better limiting behaviour for p = 2 it would be necessary to 
go to much larger values of L. Another useful check of the irrelevance of the long ranged 
boundary fields is the value of the Casimir amplitude A ++ , a universal quantity defined as 
A n (d - 1) = Wn(0) (see Eq. flTIJ-. Its value is equal to -7r/48 « 0.065 in d — 2 Ising model 



with (++) boundary conditions |42j and should be the same for all p > 2. For each value 
of p, hi w 8.16 and L = 300 we have calculated the Casimir amplitude and found a very 
good agreement with the prediction. For example, A ++ w 0.065 for p = 50, ~ 0.067 for 
p = 4, and ~ 0.072 for p = 3. Again for p = 2 corrections to the finite-size scaling become 
important and we observe the significant deviation from the universal value, i.e. A ++ ph 2.5. 
The above results show that the asymptotic behaviour of the solvation force for the Ising 

n 

system below T c agrees with the formula (0) obtained by Attard et al [30j for fluids, i.e., 
for large L, f so iv{L) decays with the same power-law as the boundary field. In order to 
enquire about the amplitude, in Fig. El we plot f so i v calculated for L = 100 and hi ~ 
8.16, scaled with its value at T/T c = 0.75 as a function of the reduced temperature T/T c , 
i.e. f* olv = f so i v (T/T c ;hi,L,p)/ f so i v (0.75;hi,L,p). It follows from this figure, that for 
T < T c the amplitude of the power-law is the same for all p and depends weakly on the 
temperature. We have checked that its value is equal to 2m*(T)hi, where m*(T) is the bulk 
spontaneous magnetization. Thus, for large L, f so i v ~ 2m*(T)hi/L p and hence f* olv (T/T c ) ~ 
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m*(T/T c )/m* (0.7 '5). Note however, m* = above T c , and in this temperature range we 
observed a different power-law decay. Moreover, the force becomes attractive which is in 
disagreement with (jSJ. In the next Section we investigate the origin of this behaviour within 
a continuum Landau theory. 

It is also instructive to observe the effect of changing the amplitude hi on the temperature 
dependence of the solvation force for various range of the surface potential. Figs.E|a,b display 
fsoiv as a function of the variable y at fixed L = 100 along the bulk coexistence line h = for 
a selection of the strength parameter hi and for p = 50 and 2, respectively. As the amplitude 
hi varies between zero and 8.16 we observe a nontrivial crossover behaviour of the solvation 
force in the critical regime, associated with the change of the position of the minimum from 
the temperature below T c to the temperature above T c . Note that for short ranged boundary 
fields hi = corresponds to the ordinary transition universality class, whereas hi ~ 8.16 
corresponds to the normal transition universality class. In the crossover between these two 
universality classes a scaling variable L/l\ where li ~ h x v Y becomes relevant. L/l\ <C 1 
corresponds to the ordinary transition and L/li ^> 1 corresponds to the normal transition. 
For L/li = 0(1) a strong deviations from the universal behaviour are expected. In d = 2 
Ising model li ~ h{ 2 . For a long ranged boundary fields there appear an additional scaling 
variable hiL( d+2 ~ r, " 2 ~ p which, as already mentioned above, is irrelevant in the RG sense but 
gives corrections to the finite size scaling which may be important for small wall separations 



4l| . In Fig. E]we can see that as hi increases from zero the minimum located below T c 



reduces and shifts towards T c . At a certain small value of hi a second shallow minimum 
appears on the other side of T c . At some small value of hi, which depends on p, the minima 
become symmetric. For p = 50, which is almost a short ranged boundary field, it takes 
place at hi ~ 0.0605, i.e. when L/l\ = 0(1). As hi increases further the minimum below 
T c diminishes and a minimum above T c becomes deeper; finally a single minimum above T c 
remains. The formation of a two maxima structure in the function f so i v (T) is also observed 
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by varying the distance between the walls at fixed small value of hi. Examination of the 
shapes of the functions of Fig. 0] reveals that the symmetry between f so i v (T) for the free and 
(++) boundary conditions given by Eq. (jl2j) is broken for the long ranged surface fields. 
This may be connected with the presence of a new scaling variable hiL( d+2 ~ v ^ 2 ~ p . Clearly 
the solvation force for hi = does not change with the parameter p, whereas for large values 
of hi the depth of the minimum of f so iv{T) increases and its position moves monotonically 
towards higher values of T as the range of the surface potential increases. 

The change of the locus of the minimum of the solvation force as the parameter p is varied 
is consistent with the behaviour of the specific heat Cn(L,T;p, hi) = —T(d 2 f/dT 2 )L and of 
longitudinal spin-spin correlation length £y. These quantities are readily calculable from the 
total free energy obtained in DMRG. £11 may be expressed in terms of the ratio of the largest 
A and second largest Ai eigenvalues of the transfer matrix (L, T; p, hi) = — ln[Ai/A ]. 
We have calculated Ch and ^ x as a function of temperature for L = 100, hi ~ 8.16 
and various p. Both quantities exhibit extrema above T c which decrease and shift towards 
higher values of T as the range of the surface potential increases (see Fig. El a, b). If the 
position of the extrema of f so iv(T), Ch(T), and (T) at h = is governed by capillary 



condensation 32[ (see Sec. |H}, then their shift towards higher values of T indicates that the 
(capillary) critical point moves further away from the bulk coexistence line h = as the 
range of the surface potential increases. 



III. ANALYSIS OF LANDAU FUNCTIONAL 

In this section we analyse the Landau theory corresponding to the Ising model considered 
in Sec. |H] in order to understand why at high temperatures the solvation force is negative 
and, for large wall separation L, decays as L~^ p+l \ whereas at low temperatures f so iv(L) is 
positive and the form of decay agrees with the prediction JSJ), i.e. f so i v (L) ~ L~ p , L — > oo. 

In the Landau theory the magnetization profile m(z) in a slit of width L subject to the 
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boundary 
per k B T) 



nc 



d H(z) is obtained by minimizing the free-energy functional (per unit area and 



F[m] 



dz 



\b + f b (m(z)) - H(z)m(z) 



(15) 



where b is a positive constant. fb(m) is the bulk free-energy density 

fb{m) = -arm 2 + -am 4 — hm, 



(16) 



where r is the reduced deviation from the bulk critical temperature, a is constant and h is 
the bulk magnetic field. We consider a long ranged boundary field H(z) of the form 



H(z) = ht 



+ 



1 



(17) 



{z + X)p (L + X-z)p_ 
where we have introduced the parameter A satisfying L ^> A ^ 1 in order to assure that 
H(z) is well behaved at the boundaries. In the final analysis we take the limit X/L — » 0. 

In order to obtain the solvation force one has to find the equilibrium magnetization profile 
m(z) eq and then calculate the total excess free energy per unit area 



f ex (L) = F(L)-Lf b {m b ), 



(18) 



where F(L) = F[m eq ] and m;, is the bulk magnetization at given T, h. The solvation force 
is defined by Eq. (JTUJ). 

Henceforward, since we are interested in the asymptotic decay of f so iv{L) for temperatures 
above T c and at vanishing bulk magnetic field, for our analysis it is sufficient to keep only the 
term quadratic in m in the expression for the bulk free energy density. In this temperature 
range the bulk correlation length = (b/ar) 1 ^ 2 and = m* = 0. 

Minimization of Eq. ()15|) yields an Euler-Lagrange equation 

d 2 m(z) 



dz 2 



Z b - 2 m(z) - (l/b)H(z) 



(19) 



We assume the following boundary conditions 



m(0) = m(L) = 1, 



(20) 
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valid for the considered case, i.e. large hi, strongly adsorbing walls. Imposing the condition 
that at equilibrium the profiles are symmetric around z = L/2 yields further 



dm 



0. 



(21) 



dz z=L/2 

The solution to the Euler-Lagrange equation (|19|) is equal to the sum of the general solution, 
rrig(z), of the corresponding homogeneous differential equation (with H(z) = 0), and any 
particular solution, m p {z), of the inhomogeneous equation (with H(z) ^ 0): 



m(z) = rrig(z) + m p (z). 



(22) 



Above T c the general solution of equation (jTUj) with H(z) = satisfying the symmetry 
condition (J2T|) is 



m g (z) = A 



+ e 



-(£-*)/& 



(23) 



The required solution of the inhomogeneous equation can be expressed in the form of a 
power series in H(z) 



m p (z) = Bi 



+ 



1 



z + xy (l + x — z)p 



+ B 2 



1 



1 



(z + X)p+ 2 (L + X-z)p+ 2 



+ ... (24) 



where the constants B\, B 2 , ■ ■ ■ have to be determined by equating coefficients. We substitute 
()24|) for m(z) and its second derivative in the differential equation ()19j) and notice that for 
large separation L — > oo the lhs of this equation is subdominant to the rhs. Therefore, in 
the limit L — > oo, we can approximate m p (z) by the solution of the equation 



= & 2 m(z)-(l/b)H(z), 



(25) 



i.e., m p (z) ~ (£%/b)H(z). Thus the general solution of (fTH|) is 



m(z) ~ A 



+ e 



+ jH{z), 



(26) 



where the constant A can be determined from the boundary conditions (|20p. Substituting 
the above solution into the functional (jl5j). performing the integral over z and then the 
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derivative with respect to L, one can find the asymptotic behaviour of f so iv{L) when L — > oo. 
Notice, that for T > T c the bulk spontaneous magnetization m* = 0, so that fb(jn b ) = in 
the Eq. ()18)) for the excess free energy per unit area. 

Substitution of the equilibrium magnetization profile ()26j) into the integrand in (fTo"]) 
yields terms purely exponentially and purely algebraically decaying with z, as well as terms 
in which an exponential and an algebraic decay mixes together. There are no terms decaying 
in the same fashion as the boundary field. Such terms, if they were to exist, would yield 
the power-law decay of f so i v that is consisitent with the result by Attard et al. (Eq. ©). 
Careful analysis reveals that the leading asymptotic behavior of the solvation force arises 
from a purely algebraically decaying term in (b/2^)m 2 (z) — H(z)m(z), namely 

e b h\ i 



b (A + z)p(L + A - z)p' 
The contribution to the f so i v from the above term is (see Appendix) 



(27) 



Wd [ L dz 1 - ^ P 1 +0(L-^) L^oc (28) 

Other purely algebraically decaying terms in {b/2^l)m?{z) — H(z)m(z) give contributions 
to the solvation force that are 0(L~( 2p+1 ^) and the mixed terms cancel out. 

Contributions to the solvation force arising from (1/2)6 (dm/dz) 2 decay faster than 
£-(p+i) f or L _> oo, i n this expression there are four mixed terms: 



(z + X)p+ v (L + X-z)p +1 ' 

e -{L-z)/£ b e ~z/ib 



(29) 
(30) 



(^ + A)p +1 ' (L + \-z)p+ 1 ' 

The asymptotic behaviour of contributions to the solvation force arising from the above 
terms is (see Appendix) 

dL Jo (L + X - z)p+i ~ dL Jo dZ (X + z)p+* ~ (L+ X) p+2 + ° {L h L °° 

(31) 
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and 

~dLJo dZ {\ + z)P+^ = ~dLJo dZ (L + X- z)v+^ = ~ (A + L)^ 1 + W$ ' L ~" °°' 

(32) 

Purely algebraically decaying terms in the expression —(b/2)m(z) (d 2 m/dz 2 ) give contribu- 
tions 0(L-^+ 2 )). 

In conclusion, we have found that for T > T c 

W£) = -§^^ + 0(L-<^>), L^oo. (33) 

Note that since b > the solvation force is negative. This asymptotic behaviour does not 
change if we take into account higher order terms in the solution for m p (z) (see Eq. (124)1 ) . 

The fact that in the limit L — > oo the solvation force decays faster than the boundary 
field is due to the absence in the expression for the free energy density of terms proportional 
to H(z). Below T c the solution of the Euler-Lagrange equation (Eq. (flUjl with an additional 
term r~ 1 ^ b ~ 2 m 3 (z) in the rhs ) for large L has the form m(z) = m* + ih(z), where m* ^ is 
the bulk spontaneous magnetization. Similarly to the case of T > T c , the function m(z) can 
be expressed as a sum of the solution of the corresponding homogeneous differential equation 
(with H(z) = 0), mh(z), and a solution of the inhomogeneous equation (with H(z) ^ 0), 
m inh(z), i.e. rh(z) = m h (z) + m inh (z) . For large L the approximation for m inh {z) is 
still valid, therefore in the expression for the free energy we can expect terms ~ m*H(z) 
which then yield an asymptotic decay of the solvation force of the same form as that of the 
boundary field. 

IV. DENSITY FUNCTIONAL THEORY RESULTS 

The model considered in this Section is a van der Waals fluid of the bulk density pb 
confined in a slit of width L. Each of the walls interacts with the fluid via the potential 
described by Eq. (J3J) with p — 2,3. The total external potential of the system V ext (z) is a 
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sum of the contributions from both walls, V ext (z) = V s (z) + V S (L — z). The fluid particles 



interact via the standard Lennard- Jones potential 

( [/ n 12 / \6 

\4e„ 

u(r) 



, r < r cut 

(34) 

, r > r cut 



where £// describes the strength of the fluid-fluid interactions and r cut is the cut-off distance. 
We set r cut = 2.5a. The system is studied by means of a density functional theory (DFT) 



44j. Within this approach the grand potential of the system is a functional of the local 
density p(r) 

Q,\p] = F[p] + J d 3 r p(r)(V ext (r) - /i) (35) 

where /i is the chemical potential of the fluid. The free energy functional F is a sum of two 
parts, F = Fid + F ex . The ideal gas contribution is known exactly 

f3F ld = J d 3 r [ln(A 3 p(r) - l]p(r) , (36) 

where (3 = (fc^T)" 1 and A is the de Broglie wavelength. The excess (over ideal) free energy 
is a sum of reference hard-sphere F^ and attractive F^ tt ' contributions. The latter is 
evaluated in a mean-field fashion 

Fjf ) = ljd 3 rj d?r> p{vW)u WCA {\v - v'\) (37) 

where uwcA( r ) corresponds to the Weeks- Chandler- Andersen |43J division of the interpar- 
ticle potential 

( \ J ~ £ff ' T<2la /OON 

UwcA{r) = < (38) 
I u{r) , r > 26(7 

The reference hard-sphere part of the excess free energy is evaluated within the framework 
of the Fundamental Measure theory (FMT) of Rosenfeld 45, 4|| 

PF<M = J d 3 r $({n a }) , (39) 

where n a denote weighted densities 

/"dVp(r') w a {r-r') , (40) 
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n n r 



with six different geometrical weight functions (four scalar and two vector-like 
There are several expressions for the excess free energy density $. We have chosen for the 
present problem the original Rosenfeld functional, where the excess free energy density is 
given by 

A / r n n /, \ i nina-ni-na n\ - 3n 2 n 2 ■ n 2 

$ KJ = -n In 1 -n 3 + + — — — . 41 

1 — 77,3 247T(1 — ri3) z 

The density profile p(r) = p(z), for planar walls, is obtained by solving the Euler-Lagrange 
equation, i.e. from 

<*r°- (42) 

The solvation force f so iv{L) (or excess pressure) can be obtained from 

U - 4 (^) (43, 

where A denotes the area of the wall and Q ex = Q + pV is the excess grand potential of 
the system. Here p is the pressure of the reservoir at the chemical potential \x and the bulk 
density pb and V is the total volume. Statistical mechanical sum rules for a confined 
fluid lead to another expression for the solvation force 

roo Q 

fsoh, = ~P- J dzp(z) — V s (z) (44) 

We have used the above equation to check the accuracy of the numerics. 

Without loss of generality we chose the LJ a as a unit of length, and introduce the 
following reduced units, f so i v = PfsoivC 3 , T* = ^— , pi = pb<J 3 . For the system in question 
the reduced critical temperature and the reduced critical density for the gas-liquid transition 
are T* = 1.319442 and p* = 0.245736, respectively. 

We start by reporting the temperature dependence of the rescaled solvation force fs 0lv , 
i.e. the solvation force divided by its value at T/T c = 0.615, for systems with p = 3, 2 
and L = 50.4<r (see Fig- EJ) - In analogy with the results presented in Sec. Ill Cl the chemical 
potential p{T) (or, equivalently the bulk density of the reservoir, pi) is fixed such that 
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the bulk fluid is slightly off coexistence on the liquid side of the bulk coexistence curve. We 
notice that for the system with the wall-fluid interactions of the finite range, i.e. for truncated 
wall-fluid potential (jHJ with p = 3 and cut-off z = 2.5a, (see the inset) the solvation force 
is extremaly small away from T c . At low temperature the system is on the oscillatory side 
of the Fisher- Widom line, so the solvation force should decay in an exponentially-damped 



oscillatory fashion i.e. f so iv(L — > oo) ~ exp(— a L) cos(aiL) |48(. Consequently, for the large 
separations considered here the solvation force well below T c becomes extremely small in 
magnitude and perishes in the numerical noise, i.e. its magnitude is smaller than 10~ 10 . 
Close to the critical region a pronounced decrease in the solvation force is found (see the 
inset) with its minimum located slightly above T c . Note that for T > T c we follow the 
critical isochore Pb(p, r) = p c and the decay is expected to be purely exponential, at least 
in the range shown in Fig. El For this value of L the magnitude is extremely small. 

For the systems with long ranged wall-fluid potentials we find that for p = 3 (Fig. EJ 
solid line) the solvation force is positive away from the critical region. Only very close to 
T c does f* olv change sign and become negative. The minimum is located slightly above 
T c . For still stronger fluid-wall potentials i.e. for p = 2 (Fig. El dashed line) the solvation 
force is positive (repulsive) through the entire range of temperatures under consideration. 
The rescaled solvation force is almost identical for both values of p at temperatures away 
from T c which is consistent with the result of Eq. (JHJ), i.e. f* olv ~ p(T) . As the systems 
move towards the critical region, f* olv begins to be different for different powers of the wall- 
fluid potentials. However the minimum of f* olv is located, similarly to previous cases, slightly 
above the critical temperature. These findings are in accordance with general considerations 
presented in the Introduction. The minimum of the solvation force is connected with the 



critical Casimir effect 3l(. For short ranged wall-fluid potentials this effect is dominant (see 
the inset to Fig. El) as the direct influence of the wall is negligible at large separations. When 
the long ranged wall-fluid potentials are introduced, the contribution from the regular part 
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of the solvation force f™f v , which is different for different powers of the wall-fluid potential, 
dominates for temperatures away from the T c while for systems in the critical region the 
singular part of the solvation force /J„jjf comes into play. The dominant decay of the singular 
part is the same for both value of p. In the present DFT approach the critical Casimir force 
is treated in mean-field, therefore jf^f is given by (jllj) with d — 4 and the mean-field value 
of the exponent v [y MF = 1/2). 

In Figs. [7HB1 we show the L-dependence of the solvation force for systems with long- 
range wall-fluid potentials for the bulk reservoir state T* = 1.0, pi = 0.6148. Again the 
temperature and chemical potential were fixed such that the system is just slightly away from 
coexistence on the liquid branch of the coexistence curve. We observe that for intermediate L 
the well-known oscillatory behaviour of f so i v characteristic, for short separations, is damped 
and changes gradually to the power-law decay enforced by the external (long ranged) wall- 
fluid potential. The asymptotic behaviour of the solvation force is demonstrated in Fig. 0d 
and Fig. IHt", where we plot the logarithm of f so i v as a function of the logarithm of the 
inverse separation (symbols) along with best straight line fits with imposed slope (solid 
lines) for the systems with p = 3, and 2, respectively. These are similar to these in Fig. 01 
presented in Sec. Ill CI In both cases we find a nice agreement between the DFT results 
and the predictions from Eq. The asymptotic form of the decay is quite visible already 
for medium separations, i.e. for L > 50. In order to investigate further the asymptotic 
behaviour of the solvation force we have performed least-square fits assuming a power-law, 
i.e. assuming f so i v (L) = a L K with a , k taken as fit parameters. For the system presented 
in Figs. EHHl best fits performed for 190 > L > 110 are 4.69x~ 2,99 and 4.62x~ L ", respectively, 
where x = L/a. The value of the constant ao predicted by Eq. (JHJ) is 4.92 and differs from 
those given by least-square fits by ~ 6%. We also note that the constant ao should be 
the same for both potentials considered here (it depends on only the density pi and the 
parameter B = 4ef w which is same for both) and indeed, to a good approximation, this is 
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the case here. 



It has been well established jl2t Il3| that if the chemical potential is fixed such that the 
fluid is on the gas side of the liquid-gas phase diagram the solvation force can (as a function 
of the separation) exhibit a jump that is a direct manifestation of capillary condensation. 
The discontinuity can appear for both short ranged and long ranged wall-fluid potentials. If 
the wall-fluid interactions are short ranged, f so i v (L) will change from small negative (weakly 
attractive) values at large L (corresponding to the 'gas' phase) to larger negative values 



for smaller L (corresponding to the 'liquid' phase). As argued in the Ref. jl3| on the basis 
of macroscopic thermodynamic arguments there is always the term — A/i(p/ — p g ) in the 
expression for f so i v for the capillary condensed 'fluid' which gives rise to the aforementionec 
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jump. On the basis of the theory of the asymptotic decay of the correlation functions [1C , 
it was argued that the solvation force should decay (asymptotically) in the same manner 
as the bulk pair correlation function. In the present case the fluid lies on the monotonic 
side of the Fisher- Widom line. A schematic plot of f so i v is shown in Fig. (for clarity the 
oscillatory part of the solvation force, relevant for small separations is omitted here). 

When the wall-fluid interactions become long ranged the solvation force should decay 
asymptotically in the same fashion as the wall-fluid potential because its contribution to the 
fsoiv will be dominant with respect to the fast-decaying (exponential) fluid-fluid contribution. 
On the other hand, the discontinuity in f so i v connected with capillary condensation must 
remain. Thus one can anticipate that the solvation force should behave in the manner 
presented in Fig.^b (again, for clarity the oscillatory part of the solvation force, relevant for 
small separations is omitted here). Namely, upon increasing the wall- wall separation L, the 
solvation force jumps from larger to smaller but still negative values. This discontinuity is 
associated with capillary condensation. Next, f so i v changes its sign (the point denoted by Z), 
attains a maximum (the point denoted by M) followed by the inflection point (denoted by 
/) and, finally, reaches the region of the asymptotic decay dictated by Eq. (j5J). The specific 
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example of the solvation force calculated by DFT for the system with p = 3, T* = 1.2 
and pi = 0.06 shown in Fig. HU1 fits well into the general behaviour described above. The 
jump associated with the capillary condensation occuring at L ~ 10.5 (c.f. Fig. ITUk) is from 
larger to smaller but still negative values of f so iv As the wall- wall separation is increased 
(c.f. Fig.llOb) f so i v changes its sign around L ~ 72, attains a maximum at L ~ 95. After an 
inflection point at L ~ 130 the solvation force changes slowly towards its asymptotic form of 
decay. It must be mentioned however that even for the largest separation studied (L = 190) 
the asymptotics (i.e. the decay with a power-law the same as the decay of the wall-fluid 
potential) could not be still reached. We think one needs to go to the separations as large 
as several hundreds to rich the asymptotic range. 

V. SUMMARY AND CONCLUSIONS 

In this paper we have performed calculations of the solvation force for an Ising film 
subject to long ranged boundary fields (J2J) with p = 50,4,3, and 2, and for a truncated 
LJ fluid confined between two planar walls that exert a 9 — p with p = 2, 3 wall-fluid 
potential (J3j). For an Ising film results have been obtained in d = 2 for states along the line 
of the bulk two-phase coexistence h = by means of the DMRG method which takes into 
account fluctuations of the order parameter. Since in two dimensions critical fluctuations are 
particularly strong, results obtained in this model serve as an ultimate test of the effects of 
fluctuations on the behaviour of the solvation force. For LJ fluid results have been obtained 
for several states on a liquid side of the bulk two-phase coexistence and on the critical 
isochore for T > T c within a nonlocal DFT which is a mean field theory. This approach 
accounts for packing effects and hence oscillations in the solvation force. 

We observe major differences in the behaviour of the solvation force between the both 
models. In the Ising system f soiv at low temperatures is positive (repulsive) and decays 
for large L in the same fashion as the boundary field, i.e., f so i v = L~ p , whereas at high 
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temperatures f so i v is negative (attractive) and the asymptotic decay is of the higher order 
than that of the boundary field, i.e., f so i v = L~~( p+1 \ In the LJ fluid system f so i v is always 
repulsive away from the critical region and decays for large L with the the same power 
law as the wall-fluid potential, which is consistent with the general result based on the 



wall-particle Ornstein-Zernike equations |30j]. As discussed within a Landau approach in the 
Sec. Illll the origin of this discrepancy is due to the specific symmetry of the Ising model with 
the spontaneous magnetization m* equal to zero above T c ; note that for a fluid p c {T) ^ 
above T c . 

Our results imply that for a confined fluid with short ranged fluid-fluid interactions 
and long ranged wall-fluid potential the solvation force, for large wall separations, can be 
expressed as a sum of a regular part f s e f v , which decays in the same fashion as the wall- 



fluid potential, and a part f s s 1 ^ arising from L-dependent singular contribution 



to the 



solv 

free energy, which is a close analog of the Casimir force in electromagnetism. The singular 
contribution to the free energy is responsible for the critical singularities at the bulk critical 
point. Thus, for ordering field h = and L — > oo 

fsoiv ~ a + a = 2PBL~ P + (d- l)L- d W aa {rL^) (45) 

where W aa (y) is a finite-size scaling function for system with two identical walls a, and 
B = 4Ef w for a wall- fluid potential of a form (JHJ). In mean-field d takes the value 4 in the 
final term of Eq. ()45|) . Such a decomposition applies also for an Ising spins system below 
T c with pB replaced by m*h\. Above T c , m* vanishes and the leading regular part of the 
solvation force becomes ~ L~(p +1 ) as already mentioned above. The scaling function W aa (y) 
is negative and vanishingly small away from the critical region, therefore, in d = 2 and 3 and 
for all values of the parameter p the positive regular part of the solvation force dominates 
away from T c . Fluctuation induced attractive Casimir force is particularly strong in two 
dimensions and dominates the behaviour of the solvation force in the critical region for all 
considered values of the parameter p - see Fig. [T] In mean field fl™£ is much weaker and 
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only when p > d = 3 the solvation force becomes weakly negative near T c . For p = 2 the 
regular part j r s e f v gives the leading decay and although the solvation force exhibits minimum 
near T c , it remains positive for all temperatures as is seen in Fig. El 

As a last remark we note that the inclusion of power-law fluid-fluid interactions would 
modify our results. This can be infered from the asymptotic integral equation results for 
the solvation force of the full LJ fluid - see Eq. l{4*|). The final term is associated with the 
_ r -n decay of the fluid-fluid potential. For LJ n = 6 so Eq. (J3J) predicts an additional 
attractive term ~ L~ 3 in the solvation force. This should then be included in Eq. (|45j). 
This contribution competes with other terms and may lead to even richer behaviour of the 
solvation force. We want to stress that in our studies we have neglected the contribution 
due to the direct wall- wall interaction potential |51|. 

After the completion of our calculations we learnt of a recent article by Pertsin and 
Grunze [52] describing the results of Monte Carlo simulations of a Lennard- Jones fluid 
confined between two planar walls. Since their results were at complete variance with ours 
and with the general predictions of Eq. (£Q) we decided to perform DFT calculations for the 
same state point as in the simulations of Ref. 52|. The results are presented in Appendix 
B. As expected, we find no evidence for the extremely long ranged solvation force reported 



in Ref. 



We have since ascertained, via correspondence with Professor A. Pertsin and 



Professor R. Evans, that the simulations reported in Ref. 5_2j were flawed and that the 
authors shall publish an erratum pointing this out. The new results, however, are in a good 
qualitative agreement with our DFT results. 
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APPENDIX: A 



In order to evaluate contributions to the solvation force arising from terms ([27]). (|29p and 
flUJ) we use Eqs. (HUJ), (Hi), and (JTHJ). 



First we integrate over z to obtain contributions to the total free energy. We find | 

dz L ^ 2 2 k+1 (2p- l)(2p-3)---(2p- 1 -2k) 

o (A + z)p(L + A-^)p~2p-lj^(p-l)(p-2)---(p-fc - l)A fc+1 (LA + A^- 1 ' 

(A.l) 

where A = (4LA + 4A + L 2 ). In the limit of L — > oo and A/L — > the above integral 
simplifies 

p-2 

£ 

fc=0 



rf£ ^ Z 2 fc+1 (2p-3)---(2p - 1 -2fc) 

o (A + z)p(L + X-z)p ^ L ^°° £L (p - l)(p - 2) • • • (p - jfe - ^Ap-^ 1 ^' 



(A.2) 



It follows that the contribution to the solvation force arising from the term (|27|) is 



dz- 



1 



6 dLJo (X + z)p(L + X - z)p 
Integration of terms (|29|) yields 



L e -{L-z)/£ b 



dz- 



and 



o (L + A - z)p Jo (A - z) 



i p- 1 



(A.4) 



(£ + A-s)* " ^(p-l)!^ 

1 



A fc (L + A) fc 



+ 



?6 £b 



(A.5) 



where Ei(x) is the Exponential- Integral Function. 

An asymptotic representation of the function Ei(x) for large x is 



x — > oo 



k=l 



(A.6) 
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Thus, in a limit of large L, the second term in the expression for the integral (jA.5|) behaves 



as 



1 
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„ . / A s _ , A . 
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+ 0((L + A)-<" +1) ) (L — * oo). (A. 7) 



Using the above asymptotic expression we obtain the asymptotic behaviour of the integral 
31) in a limit of L — > oo 
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It follows that 
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Now, consider contributions to the solvation force arising from terms 



We have 
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APPENDIX: B 



In order to compare with the results of Ref. 52| we carriec 
T* = 1.0 and pi = 0.6853, the simulation state point. In Ref. (^2] 



out DFT calculations for 
the LJ potential is cut-off 
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at r cut = 2.5a, the same as in our DFT, and the wall-fluid potentials have the same form 
as ©• Fig. ITTk shows our results for f so i v {L) for the case of a 9-3 wall-fluid potential, i.e. 
p = 3, whilst Fig. ITTb refers to the truncated 9-4 wall-potential, set to zero at z — 2.5a. In 
the first case we confirmed that for large L, f so iv{L) is positive and decays as L~ 3 . Indeed the 
fit to this power law is as good as that for the neighbouring state point shown in Fig. [7| For 
the truncated 9-4 potential f so iv(L) exhibits the expected exponentially damped oscillatory 
decay for large L. Thus our results for this state point are also in agreement with the general 
theory described in the Introduction. By contrast Pertsin and Grunze 52[ find for large L 
an attractive solvation force which decays very slowly, roughly as L _1 , for both the 9-3 and 
the truncated 10-4 wall- fluid potentials (see their Fig. 3a). Such behaviour would seem to 
be unphysical. 
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Figure Captions 

Fig.l Solvation force (in units of J) for an Ising strip of the width L = 300 at vanishing 
bulk magnetic field subject to long ranged boundary fields Hi = Hf + Hl +1 _ b Hf = 
hi/l p as a function of the scaling variable y = r(V) x l v = tL for various values of p 
and the amplitude hi ~ 8.16 calculated using DMRG method. The inset shows, on 
an expanded scale, that results for p = 50 and p = 4 differ only very slightly in the 
near-critical region. For all values of p f so i v is attractive above T c and, except for 
p = 50, repulsive below T c . 

Fig.2 Log-log plot of the solvation force (in units of J) versus 1/L for two values of 
the power p describing the decay of the boundary field iff, p = 2 and 4. (a) displays 
DMRG results for T = 0.79T C (filled symbols) and T = T c (unfilled symbols). Straight 
lines correspond to f so i v (L) ~ L~ p . Data at T c align parallel to the line of a slope 
2 showing the result f so i v {L) ~ (see Eq. (1)) with d = 2. (b) displays DMRG 
results for T = 1.23T C . Straight lines correspond to f so iv{L) ~ L~^ p+1 \ 

Fig.3 Solvation force for the Ising strip of width L = 100, hi ~ 8.16 rescaled with its 
value at T/T c = 0.75, plotted as a function of the reduced temperature. Below T c , 
the results calculated for three different values of the power p lie on a common curve 
~ m*(T/T c ), where m*(T) is bulk spontaneous magnetization. This behaviour agrees 
with the prediction Eq. (JHJ) with 2pB replaced by 2m*(T)hi. 

Fig.4 Solvation force (in units of J) for the same system as in Fig. El as a function 
of the variable y for a selection of values of the strength parameter hi and for two 
values of the parameter p: (a) p = 50 and (b) p — 2. Results for different values 
of hi in Fig. 0Jd are represented by the same symbols as in Fig. H^,. It is seen that 
the symmetry, Eq. (fT2"|) between f so i v (y) for the free (hi = 0) and (++) (hi ~ 8.16) 
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boundary conditions breaks when the boundary fields become long ranged. For weak 
hi the solvation force exhibit two minima. 

Fig.5 (a) Specific heat Ch (in units of ks) and (b) the inverse longitudinal spin-spin 
correlation length for the same system as in Fig. El as a function of the scaling 
variable y for h% ~ 8.16 and various values of the parameter p (represented by the 
same symbols for both quantities). 

Fig. 6 Solvation force rescaled with its value at T/T c = 0.615, f* olv = 
fsoiv(T/T c )/ f so i v (0.615), as a function of the temperature for systems with long ranged 
fluid-wall potentials for L = 50.4a evaluated from density functional theory. The 
chemical potential fi(T) is chosen so that the reservoir density bulk density pl(T) is 
close to that of the bulk coexisting liquid for T < T c and pb(p,,T) = p c for T > T c . 
Solid line corresponds to the system with p = 3, while dashed line is for the system 
with p = 2. The inset shows the temperature dependence of f so i v for the system with 
fluid-wall interactions of the finite range. 

Fig. 7 The oscillatory (part a) and asymptotic shown on a log — log plot (part b) regions 
of the solvation force for the system with p = 3 evaluated from density functional 
theory. The temperature (T* = 1.0) and the reservoir density (pi = 0.6148) are 
fixed such that the fluid is just slightly off the coexistence on the liquid branch of the 
coexistence curve. 

Fig. 8 The oscillatory (part a) and asymptotic shown on a log — log plot (part b) regions 
of the solvation force for the system with p = 2 evaluated from density functional 
theory. The temperature (T* = 1.0) and the reservoir density (pi = 0.6148) are 
fixed such that the fluid is just slightly off the coexistence on the liquid branch of 
the coexistence curve. Note that for p = 2 the solvation force has more pronounced 
oscillations than for p = 3. 
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Fig. 9 Schematic plot of the solvation force for systems with short ranged (part a) and 
long ranged (part b) fluid-wall interactions in the presence of capillary condensation. 
Dots denoted by Z, M and I in part b mark characteristic points of the solvation force: 
the zero, the maximum and the inflection point, respectively. 

Fig. 10 The short distance (part a) and asymptotic (part b) regions of the solvation 
force for the system with long ranged (p = 3) wall-fluid interactions. Note the change 
of scale in each part. Capillary condensation gives rise to the jump at L — 10.5. 
For L < 10.5 the confined fluid is a 'liquid' which can exhibit oscillations. The bulk 
reservoir corresponds to (T* = 1.2, pi = 0.06). 

Fig. 11 Solvation force for the LJ fluid with the short ranged and long ranged wall-fluid 
potentials calculated within DFT for T* = 1.0 and pi = 0.6853, the same state point 



as in the simulations by Pertsin and Grunze 5_2j. (a) is for the full 9-3 LJ potential. 
Oscillations characteristic for short wall-wall separations are gradually damped and 
changed into the power law asymptotic decay, (b) is for the truncated 9-4 LJ potential 
with the cut-off distance 2.5a. Oscillations around are present for all wall-wall 
separations. 
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FIG. 1: Solvation force (in units of J) for an Ising strip of the width L = 300 at vanishing bulk 
magnetic field subject to long ranged boundary fields Hi = Hf + H S L+1 _^ Hf = h\/l p as a function 
of the scaling variable y = t(V) x I v = tL for various values of p and the amplitude h\ w 8.16 
calculated using DMRG method. The inset shows, on an expanded scale, that results for p = 50 
and p = 4 differ only very slightly in the near-critical region. For all values of p f so i v is attractive 
above T c and, except for p = 50, repulsive below T c . 
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FIG. 2: Log-log plot of the solvation force (in units of J) versus 1/L for two values of the 
power p describing the decay of the boundary field Hf, p = 2 and 4. (a) displays DMRG results 
for T = 0.79T C (filled symbols) and T = T c (unfilled symbols). Straight lines correspond to 
f so l v (L) ~ L~ p . Data at T c align parallel to the line of a slope 2 showing the result f so iv(L) ~ L~ d 
(see Eq. (1)) with d = 2. (b) displays DMRG results for T = 1.23T C . Straight lines correspond to 
f solv (L) ~ L-(f +1 ). 
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FIG. 3: Solvation force for the Ising strip of width L = 100, h\ ~ 8.16 rescaled with its value at 
T/T c = 0.75, plotted as a function of the reduced temperature. Below T c , the results calculated 
for three different values of the power p lie on a common curve ~ m*(T/T c ), where m*(T) is bulk 
spontaneous magnetization. This behaviour agrees with the prediction Eq. Q with 2pB replaced 
by 2m*(T)hx. 
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FIG. 4: Solvation force (in units of J) for the same system as in Fig.|2]as a function of the variable 
y for a selection of values of the strength parameter hi and for two values of the parameter p: (a) 
p = 50 and (b) p = 2. Results for different values of hi in Fig. are represented by the same 
symbols as in Fig. 0^. It is seen that the symmetry, Eq. (fT2*)) between f so iv(v) for the free (hi = 0) 
and (++) (hi ~ 8.16) boundary conditions breaks when the boundary fields become long ranged. 
For weak hi the solvation force exhibit two minima. 
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FIG. 5: (a) Specific heat Cjj (in units of ks) and (b) the inverse longitudinal spin-spin correlation 
length ^y 1 for the same system as in Fig. [31 as a function of the scaling variable y for hi ~ 8.16 
and various values of the parameter p (represented by the same symbols for both quantities). 
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FIG. 6: Solvation force rescaled with its value at T/T c = 0.615, f* olv = f S olv(T/T c )/f so i v (0.615), as 
a function of the temperature for systems with long ranged fluid-wall potentials for L = 50.4<r eval- 
uated from density functional theory. The chemical potential n{T) is chosen so that the reservoir 
density bulk density pl(T) is close to that of the bulk coexisting liquid for T < T c and pb{p, T) = p c 
for T > T c . Solid line corresponds to the system with p = 3, while dashed line is for the system 
with p = 2. The inset shows the temperature dependence of f so i v for the system with fluid-wall 
interactions of the finite range. 
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FIG. 7: The oscillatory (part a) and asymptotic shown on a log — log plot (part b) regions of 
the solvation force for the system with p = 3 evaluated from density functional theory. The 
temperature (T* = 1.0) and the reservoir density (pi = 0.6148) are fixed such that the fluid is just 
slightly off the coexistence on the liquid branch of the coexistence curve. 
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FIG. 8: The oscillatory (part a) and asymptotic shown on a log — log plot (part b) regions of 
the solvation force for the system with p = 2 evaluated from density functional theory. The 
temperature (T* = 1.0) and the reservoir density (pi = 0.6148) are fixed such that the fluid is just 
slightly off the coexistence on the liquid branch of the coexistence curve. Note that for p = 2 the 
solvation force has more pronounced oscillations than for p = 3. 
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FIG. 9: Schematic plot of the solvation force for systems with short ranged (part a) and long 
ranged (part b) fluid-wall interactions in the presence of capillary condensation. Dots denoted by 
Z, M and / in part b mark characteristic points of the solvation force: the zero, the maximum and 
the inflection point, respectively. 
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FIG. 10: The short distance (part a) and asymptotic (part b) regions of the solvation force for 
the system with long ranged (p = 3) wall-fluid interactions. Note the change of scale in each part. 
Capillary condensation gives rise to the jump at L = 10.5. For L < 10.5 the confined fluid is a 
'liquid' which can exhibit oscillations. The bulk reservoir corresponds to (T* = 1.2, pi = 0.06). 
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FIG. 11: Solvation force for the LJ fluid with the short ranged and long ranged wall- fluid potentials 
calculated within DFT for T* = 1.0 and = 0.6853, the same state point as in the simulations by 



Pertsin and Grunze 



521 ] . (a) is for the full 9-3 LJ potential. Oscillations characteristic for short 



wall-wall separations are gradually damped and changed into the power law asymptotic decay, 
(b) is for the truncated 9-4 LJ potential with the cut-off distance 2.5a. Oscillations around are 
present for all wall- wall separations. 
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